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Abstract 

A  one-dimensional,  non-isothermal,  two-phase  transient  model  has  been  developed  to  study  the  transient  behaviour  of  water  transport  in  the 
cathode  gas  diffusion  layer  of  PEM  fuel  cells.  The  effects  of  four  parameters,  namely  the  liquid  water  saturation  at  the  interface  of  the  gas  diffusion 
layer  and  flow  channels,  the  proportion  of  liquid  water  to  all  of  the  water  at  the  interface  of  the  cathode  catalyst  layer  and  the  gas  diffusion  layer, 
the  current  density,  and  the  contact  or  wetting  angle,  on  the  transient  distribution  of  liquid  water  saturation  in  the  cathode  gas  diffusion  layer  are 
investigated.  Especially,  the  time  needed  for  liquid  water  saturation  to  reach  steady  state  and  the  liquid  water  saturation  at  the  interface  of  the 
cathode  catalyst  layer  and  gas  diffusion  layer  are  plotted  as  functions  of  the  above  four  parameters.  The  ranges  of  water  vapour  condensation  and 
liquid  water  evaporation  are  identified  across  the  thickness  of  the  gas  diffusion  layer.  In  addition,  the  effects  of  the  above  four  parameters  on  the 
steady  state  distributions  of  gas  phase  pressure,  water  vapour  concentration,  oxygen  concentration  and  temperature  are  also  presented.  It  is  found 
that  increasing  any  one  of  the  first  three  parameters  wifi  increase  the  water  saturation  at  the  interface  of  the  catalyst  layer  and  gas  diffusion  layer,  but 
decrease  the  time  needed  for  the  liquid  water  saturation  to  reach  steady  state.  When  the  liquid  water  saturation  at  the  interface  of  the  gas  diffusion 
layer  and  flow  channels  is  high  enough  (>0.1),  the  liquid  water  saturation  at  steady  state  is  almost  uniformly  distributed  across  the  thickness  of 
the  gas  diffusion  layer.  It  is  also  found  that,  under  the  given  initial  and  boundary  conditions  in  this  paper,  evaporation  takes  place  within  the  gas 
diffusion  layer  close  to  the  channel  side  and  is  the  major  process  for  water  phase  change  at  low  current  density  (<2000 Am-2);  condensation 
occurs  close  to  the  catalyst  layer  side  within  the  gas  diffusion  layer  and  dominates  the  phase  change  at  high  current  density  (>5000  Am-2).  For 
hydrophilic  gas  diffusion  layers,  both  the  time  needed  for  liquid  water  saturation  to  reach  steady  state  and  the  water  saturation  at  the  interface  of 
the  catalyst  layer  and  gas  diffusion  layer  will  increase  when  the  contact  angle  increases;  but  for  hydrophobic  gas  diffusion  layers,  both  of  them 
decrease  when  the  contact  angle  increases. 

©  2005  Elsevier  B.Y.  All  rights  reserved. 
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1.  Introduction 

Polymer  electrolyte  membrane  (PEM)  fuel  cells  have  been 
demonstrated  as  a  promising  alternative  power  source  for 
portable,  automotive,  and  stationary  applications.  In  order  to 
commercialize  PEM  fuel  cells  several  technological  problems 
have  to  be  overcome.  Water  management  is  one  of  the  most  crit¬ 
ical  and  challenging  problems  because  insufficient  amount  of 
liquid  water  reduces  significantly  proton  conductivity  and  too 
much  water  hinders  the  reactants  coming  through  the  gas  diffu¬ 
sion  layer  (GDL)  to  the  reaction  sites  in  the  catalyst  layer  (CL) 
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causing  mass  transport  limitation.  A  fine  balance  of  liquid  water 
can  only  be  established  based  on  a  better  understanding  of  water 
transport  phenomena. 

In  the  past  decade,  two-phase  mathematical  models  for  PEM 
fuel  cells  or  their  cathodes  have  been  published  [1-15].  He 
et  al.  [1]  proposed  a  two-dimensional,  two-phase,  steady-state, 
isothermal  model  for  evaluating  the  effects  of  liquid  water  trans¬ 
port  on  the  performance  of  the  cathodes  of  PEM  fuel  cells 
with  interdigitated  flow  fields.  An  interfacial  mass  transfer  rate 
was  introduced  in  their  model  to  describe  the  mass  change 
between  the  gas  phase  and  liquid  phase  of  water.  Later,  Lin  et 
al.  [2]  further  refined  above  model  by  incorporating  a  thin  film- 
agglomerate  model  for  catalyst  layers.  Janssen  [3]  presented  a 
steady  state,  two-phase  and  two-dimensional  model  for  PEM 
fuel  cells  based  on  the  concentrated  solution  theory.  Wang  et 
al.  [4]  developed  a  two-dimensional,  isothermal  model  for  the 
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Nomenclature 

Ci  specific  heat  capacity  of  phase  i  (J  kg-1  K-1) 

Cj  mass  fraction  of  species  j  in  phase  i 

Dj  bulk  diffusion  coefficient  of  species  j  in  phase  i 
(m2  s-1) 

Dl  eff  effective  diffusion  coefficient  of  species  j  in  phase 
i  (m2  s-1) 

F  Faraday’s  constant  (C  mol-1) 

Hy ap  enthalpy  of  vaporization  of  water  (J  kg-1) 

ic  current  density  (A  m-2) 

J(s)  Leverett  function 

kCond  water  vapour  condensation  rate  (s- 1 ) 
ki  thermal  conductivity  in  phase  i  (W  m-1  K-1) 
k-  thermal  conductivity  of  species  j  in  phase  i 
(Wm-1K-1) 

kri  relative  permeability  of  phase  i 

kvap  liquid  water  evaporation  rate  (N-1  m2  s-1) 

K  permeability  of  the  GDL  (m2) 

L  thickness  of  the  cathode  GDL  (m) 

mH2°  interfacial  mass  transfer  rate  of  water 
(kg  m-3  s-1) 

Mg  average  molar  mass  of  the  gas  phase  (kg  mol-1) 
M7  molar  mass  of  species  j  (kg  mol-1 ) 

pc  capillary  pressure  (N  m-2) 

Pi  pressure  in  phase  i  (N  m-2) 

R  gas  constant  (J  mol- 1  K- 1 ) 

RH  relative  humidity 

.v  liquid  water  saturation 

t  time  coordinate  (s) 

/steady  time  needed  for  water  saturation  to  reach  steady 
state 

T  temperature  (K) 

Ui  velocity  in  phase  /(ms-1) 

v  position  coordinate  (m) 

xj  molar  fraction  of  species  j  in  phase  / 

xsat°  saturation  molar  fraction  of  water  vapour 

Greek  symbols 

a  net  drag  coefficient  in  the  membrane 

ft  proportion  factor 

A  S^{  entropy  production  for  the  cathode  (J  mol  - 1  K- 1 ) 

At  time  step  used  in  the  numerical  simulation 
£  porosity  of  the  cathode  GDL 

fi,  £2  parameters  in  current  density  expression  (Eq. 

(42))  (Am-2  V-1) 

r/c  overpotential  in  the  cathode  catalyst  layer 

0C  contact  or  wetting  angle 

Pi  dynamic  viscosity  of  phase  /  (kg  m-1  s-1) 

pi  density  of  phase  /  (kg  m-3) 

o  surface  tension  (N  m- 1 ) 

crc  electric  conductivity  of  solid  phase  (S  m-1) 


Subscripts 

c  cathode,  capillary,  or  contact 

cond  condensation 

eff  effective  value  of  a  parameter 

g  gas  phase 

j  species  index 

1  liquid  phase 

L  at  x  =  L 

r  relative 

s  solid  phase 

sat  saturation 

steady  at  steady  state 

vap  vaporization 

0  at  v  =  0  or  t  =  0 

Superscripts 
H2O  water 

i  phase  index 

N2  nitrogen 

O2  oxygen 

Pt  platinum 


air  cathode  of  PEM  fuel  cells  based  on  the  multiphase  mixture 
formulation  developed  by  Wang  and  Cheng  [5].  The  single-  and 
two-phase  regimes  for  the  water  distribution  and  transport  are 
classified  by  the  threshold  current  density  corresponding  to  the 
first  appearance  of  liquid  water  at  the  membrane/cathode  inter¬ 
face.  A  similar  model  was  published  by  You  and  Liu  [6]  to  study 
the  effects  of  operating  parameters  on  the  two-phase  mass  trans¬ 
port  in  the  cathode  of  PEM  fuel  cells.  Mazumder  and  Cole  [7] 
proposed  a  three-dimensional,  isothermal  mathematical  model 
for  PEM  fuel  cells  to  predict  the  liquid  water  transport  where  the 
water  phase  change  was  expressed  by  the  same  way  as  in  [1,2]. 
Berning  and  Djilali  [8]  also  developed  a  three-dimensional,  two- 
phase  model  for  the  cathode  and  anode  of  PEM  fuel  cells,  and 
overall  water  condensation  in  both  cathode  and  anode  was  pre¬ 
dicted. 

All  of  above  mentioned  two-phase  models  assumed  that 
the  GDL  is  completely  hydrophilic  (zero  contact  angle).  In 
addition,  they  are  mainly  steady  state  and  isothermal  mod¬ 
els.  The  liquid  water  transport  in  a  hydrophobic  GDL  was 
explained  by  Nam  and  Kaviany  in  their  one-dimensional, 
two-phase  model  [9].  Recently,  Pasaogullari  and  Wang  [10] 
extended  the  model  proposed  in  [4]  to  account  for  both  wet¬ 
ting  properties,  hydrophilicity  and  hydrophobicity,  of  GDLs. 
Weber  et  al.  [11]  discussed  the  effect  of  the  hydrophilicity  or 
hydrophobicity  of  GDLs  on  the  cell  performance  through  a  one¬ 
dimensional,  two-phase,  steady-state  model.  The  same  issue  was 
also  addressed  in  the  work  published  by  Pasaogullari  and  Wang 
[12]. 

Water  phase  change  in  PEM  fuel  cells  depends  on  pressure 
and  temperature  as  well.  A  temperature  difference  ~6°C  at 
1  A  cm-2  along  the  MEA  thickness  was  predicted  and  validated 
by  the  non-isothermal  and  two-phase  PEM  fuel  cell  model  devel- 
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oped  by  Noponen  et  al.  [13],  in  which  the  contact  resistance  for 
both  heat  transfer  and  electric  current  were  considered.  Based 
on  this  model,  Birgersson  et  al.  [14]  compared  the  simulation 
results  between  the  two  different  expressions  of  the  capillary 
pressure,  one  was  given  by  Wang  et  al.  [4]  and  the  other  was 
defined  by  Natarajan  and  Nguyen  [15]. 

The  transient  analysis  of  water  transport  through  fuel  cells 
is  needed  when  we  deal  with  the  rapid  variation  of  loads  [16], 
which  is  relevant  to  automotive  applications.  Unfortunately,  only 
few  open  reports  [16-19]  can  be  found  on  this  issue.  Ceraolo  et 
al.  [16]  investigated  the  static  and  dynamic  behaviour  of  cell  volt¬ 
age  of  PEM  fuel  cells  through  a  one-dimensional,  single  phase, 
isothermal  model.  Amphlett  et  al.  [17]  gave  an  analytical  model 
to  predict  the  transient  responses  of  a  PEM  fuel  cell  stack  dur¬ 
ing  start-up,  load  step-up  and  shut-down.  Natarajan  and  Nguyen 
[15]  developed  a  two-dimensional,  two-phase,  transient  model 
for  the  cathode  of  PEM  fuel  cells  with  conventional  gas  distrib¬ 
utors  and  the  transient  profiles  of  oxygen  concentration,  local 
current  density,  and  liquid  water  saturation  were  presented.  But 
their  model  did  not  consider  the  temperature  change  in  the  cell 
and  therefore  was  isothermal.  This  model  was  later  expanded  to  a 
pseudo-three-dimensional  one  [18]  to  account  for  the  dimension 
along  with  the  channel  direction  but  only  steady  state  analysis 
was  conducted  there.  Recently,  Wang  and  Wang  [19]  published  a 
three-dimensional,  isothermal  transient  model  to  study  the  tran¬ 
sient  dynamics  of  PEM  fuel  cell  operation  and  concluded  that 
the  time  for  fuel  cell  to  reach  steady  state  is  in  the  order  of  10  s 
due  to  the  effect  of  water  accumulation  in  the  membrane.  Gas 
transport  reaches  steady  state  after  0.1  s.  Unfortunately  it  is  a 
single-phase,  isothermal  model  and  cannot  predict  liquid  water 
transport. 

The  purpose  of  this  paper  is  to  investigate  the  transient 
behaviours  of  liquid  water  transport  in  the  cathode  GDL  of 
PEM  fuel  cells  by  means  of  a  one-dimensional,  two-phase,  non- 
isothermal  model.  It  can  determine  where  the  phase  change 
takes  place  in  the  cathode  GDL  and  under  what  conditions. 
A  parameter  is  adopted  to  account  for  the  proportion  of  liq¬ 
uid  water  to  all  of  the  water  introduced  by  electro-osmotic  drag 
and  the  electrochemical  reaction  in  the  cathode  CL  at  the  cath¬ 
ode  CL/GDL  interface.  The  effects  of  this  parameter,  as  well  as 
other  three  parameters  (the  liquid  water  saturation  at  the  interface 
of  GDL/channels,  the  working  current  density,  and  the  wetting 
or  contact  angle)  on  the  transient  and  steady  state  behaviour  of 
liquid  water  are  presented  and  discussed. 

2.  Mathematical  model 

2.7.  Basic  equations 

The  classical  multiphase  approach  developed  by  Abriola  and 
Pinder  [20]  provides  a  full  system  of  governing  equations  for 
velocities,  scalar  pressures,  scalar  liquid  saturations,  mass  con¬ 
centrations,  and  temperature.  The  governing  equations  were 
summarized  in  [5]. 

A  schematic  of  the  cathode  gas  diffusion  layer  (GDL)  of  a 
PEM  fuel  cell  is  given  in  Fig.  1.  Applying  above  multiphase 
model  to  the  cathode  GDL  where  three  species  (water,  oxygen 


MEM  CL  GDL  CL 


O2,  N2 

H20 


Lig.  1.  The  schematic  of  a  cathode  GDL  of  a  PEM  fuel  cell. 

and  nitrogen)  and  two  phases  (liquid  and  gas)  exist,  the  govern¬ 
ing  equations  can  be  obtained  as  follows. 

(i)  Mass  conservation  of  phases 


The  mass  balance  in  each  phase  is  considered  separately. 
In  the  liquid  phase 


9(p  m) 

dx 


where  s  is  the  GDL  porosity,  p\  is  the  liquid  water  density,  s  is 
the  liquid  water  saturation,  u\  is  the  liquid  water  velocity,  mH2° 
is  the  interfacial  mass  transfer  rate  of  water  from  vapour  phase 
to  liquid  phase,  and  t  and  v  are  the  time  and  position  coordinates, 
respectively. 

In  the  gas  phase 


JW  ,  d(pgUg)  _  ^  h20 

c  "  "  —  tn 


dt 


dx 


here  pg  and  ug  are  the  gas  density  and  velocity,  respectively. 


(ii)  Momentum  conservation 


The  momentum  conservations  for  the  liquid  phase  and  the  gas 
phase  are  described  by  Darcy’s  law. 

In  the  liquid  phase 


Kk ri  dp\ 

U\  = - 

fJL  1  dx 

In  the  gas  phase 

Kkrp  dp? 

Ug  = - -  - - 

fig  dx 


(3) 

(4) 


where  K  is  the  absolute  permeability  of  the  GDL,  kY\  and  kYg 
are  the  relative  permeabilities  of  liquid  phase  and  gas  phase, 
respectively,  fi\  and  fig  are  the  dynamic  viscosities  of  liquid 
phase  and  gas  phase,  respectively,  and p\  and pg  denote  the  partial 
pressure  of  liquid  phase  and  gas  phase,  respectively. 

The  difference  between  the  gas  phase  pressure  and  liquid 
phase  pressure  is  known  as  capillary  pressure, 


Pc  =  Pg~  P\ 
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(iii)  Mass  conservation  of  species 
Water  in  the  liquid  phase 

=  (6) 

Water  in  the  gas  phase 


Sjt(pgil  -  s)Cf°)  +  lx(pgugcf°) 


where  C^2°  and  C^2°  stand  for  the  mass  concentration  of  water 
in  the  liquid  phase  and  in  the  gas  phase,  respectively,  and  D^2° 
and  D^2°  are  the  diffusion  coefficients  of  water  in  the  liquid 
phase  and  in  the  gas  phase  respectively.  The  term,  7H2°,  denotes 
the  interphase  species  transfer  rate  caused  by  chemical  nonequi¬ 
librium  and/or  phase  change  at  the  interfaces  between  the  liquid 
phase  and  gas  phase. 


Oxygen 

Because  there  is  no  phase  change,  the  mass  conservation 
equation  for  oxygen  is  given  only  in  the  gas  phase 


where  C®2  is  the  mass  concentration  of  oxygen  in  the  gas  phase, 
and  D® 2  is  the  diffusion  coefficient  of  oxygen  in  the  gas  phase. 


(iv)  Energy  conservation 


respectively,  and  q\,  qg  and  qs  are  the  interphase  heat  fluxes  asso¬ 
ciated  with  liquid  phase,  gas  phase  and  solid  phase,  respectively. 

The  phase  enthalpies,  h\,  hg  and  hs ,  are  related  to  the  common 
temperature  T  via 


h\  =  f  ci  dT  +  h® 

Jo 

hg  =  J  cg  dT  +  h® 


(12) 

(13) 


T 

hs=  I  csdT  +  h°s  (14) 

Jo 


where  ci,  cg  and  cs  represent  the  effective  specific  heat  capaci¬ 
ties  of  liquid  phase,  gas  phase  and  solid  phase,  respectively,  h®, 
h®  and  h®  are  the  enthalpies  at  absolute  zero  in  three  phases, 
respectively. 


2.2.  Transformed  equations  with  initial  and  boundary 
conditions 


The  above  basic  conservation  laws  provide  a  full  system  of 
governing  equations  for  the  unknown  velocities  ( u\  and  ug),  pres¬ 
sures  (p\  and  pg ),  the  liquid  saturation  (s),  mass  concentrations 
(C^2°  and  Cg2),  and  the  common  temperature  (7). 

In  order  to  solve  the  unknowns  effectively,  the  above  govern¬ 
ing  equations  are  simplified  and  transformed  into  the  following 
forms. 


2.2.7.  Liquid  water  transport 

Substituting  the  momentum  conservation  equation  for  liquid 
water,  Eq.  (3),  into  Eq.  (1)  and  noting  that  the  partial  pressure 
in  the  liquid  phase  can  be  expressed  into  the  difference  between 
the  gas  phase  pressure  and  capillary  pressure  through  Eq.  (5), 
the  governing  equation  for  the  liquid  water  saturation  can  be 
rewritten  as 


In  GDLs,  energy  is  transferred  in  three  phases,  liquid  phase,  gas 
phase  and  solid  phase.  The  conservation  equation  in  each  phase 
is  given  as  follows. 

In  the  liquid  phase 

d  d  d  f  dT 

s-(pishi)  +  —(pmhi)  =  —  [ski  — 
at  ox  dx  \  dx 


In  the  gas  phase 


d  d 

-  s)he)  +  — (PgMg/ig) 
dt  dx 


d 

dx 


+  <?g 

(10) 


In  the  solid  phase 


(1 


d 

0—(Pshs ) 
dt 


+  qs 


(ID 


where  local  thermal  equilibrium  among  phases  has  been 
assumed  (T\  =  Tg  =  Ts  =  T),  k\,  kg  and  ks  represent  the  effective 
thermal  conductivities  of  liquid  phase,  gas  phase  and  solid  phase, 


9s  9 
sp\  —  T  — 
dt  dx 


pi 


Kkv\  d  pc  ds  Kkri  dpg 


Pi 


h2o 

=  m 


fi\  ds  dx  '  "  fi\  dx 
with  the  following  initial  and  boundary  conditions 


(15) 


s\t= 0  =  SO 


(16) 


KkY i  d  pc  ds 
pil  ds  dx 


Kkr\  dpg 

—  pi - - 

pi  dx 


x=0 


S-(2  +  4a)MH2° 

(17) 


S\X=L  =  SL  (18) 

where  so  and  sl  are  given  liquid  water  saturations,  ic  is  the  current 
density,  F  is  the  Faraday’s  constant,  MH2°  is  the  molar  mass  of 
water,  a  is  the  net  drag  coefficient,  including  both  the  electro- 
osmotic  drag  and  back  diffusion  in  the  membrane. 

A  parameter,  denoted  by  (0  <  <  1),  is  employed  in  Eq. 

(17)  to  account  for  the  liquid  water  fraction  to  the  total  amount 
of  water  (including  water  vapour)  at  the  interface  of  cathode 
CL/GDL,  which  is  introduced  by  the  electro-osmotic  drag  and 
the  electrochemical  reaction  in  the  cathode  catalyst  layer. 
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2.2.2.  Gas  transport 

The  gas  phase  transport  is  characterized  by  its  partial  pres¬ 
sure  and  the  governing  equation  can  be  obtained  as  following: 
Assuming  that  the  gas  phase  obeys  the  ideal  gas  law,  so  the  gas 
density  is  given  by 

PgMg 


Pg  — 


RT 


(19) 


where  Mg  is  the  molar  mass  of  the  gas  mixture,  R  is  the  gas 
constant  and  T  is  the  common  temperature. 

Substituted  from  Eqs.  (1),  (4)  and  (19),  Eq.  (2)  can  be  con¬ 
verted  into  a  differential  equation  of  the  gas  phase  pressure 


Mg  dpo  d 
e(l  -s)— + 


RT  dt 


Kk rg  Pg  Mg  dpg 


+ 


1 


du\ 

pi - m 

pi  [  9x 


dx 
h2o 


Jig 

Mo 


RT 


RT 


Pg  = 


—m 


dx 
h2o 


(20) 


with  the  following  initial  and  boundary  conditions 
Pg\t=0  —  PgX 


(21) 


KkVg  Pg  Mg  dpg 


_  tig 


RT  dx 


x=0 


=  — —  M°2  +  (1  -  P)  — (2  +  4a)MH2° 
4  F  P  4  F 


Pg \x=L  ~  Pg,L 

here  M°2  is  the  molar  mass  of  oxygen. 


(22) 


(23) 


2.2.3.  Water  vapour  transport 

The  liquid  water  density  is  assumed  constant  because  its 
variation  in  the  temperature  interval  of  300-363  K  is  less  than 
0(  10-2)  [14].  The  concentration  or  mass  fraction  for  liquid 
water  is  1.  So  both  the  spatial  and  temporal  derivatives  of  liquid 
water  concentration  are  0.  Eq.  (6)  is  reduced  to 

h2och2o  =  yH2o  (24) 


m 


Using  Eq.  (7),  the  governing  equation  for  water  vapour  transport 
is  derived  as 

3C»2°  9 

£pg(l— s) — r - h 


dt 


dx 


-spg(\-s)D^—f—+pgugCf° 


(/TgUg)  mn2o 
dx 


C/°  = 


(25) 


with  the  following  initial  and  boundary  conditions 
Cg2°|/=o  =  Cf0 


o^H20 

■epg(l  -  s)D^ff^—  +  pgugC^2° 


(26) 


(1  -  p)S-(2  +  4a)MH2° 


x=0 


c?2°L=l  = 


where  is  the  effective  diffusion  coefficient  of  water 

vapour. 

2.2.4.  Oxygen  transport 

Due  to  the  low  solubility  of  oxygen  in  liquid  water,  the  oxygen 
concentration  in  liquid  water  is  neglected  [6].  The  governing 
equation  for  oxygen  transport  is, 


dC°2 

spg(l  -  j)-tt-  + 
ot 


d 

dx 


dC°2 

-spg(l  -  s) Dg2eff  f  +  PgMgC°2 


,  mH20 
dx 


C ?2  =  0 


(29) 


here  D^2ff  is  the  effective  diffusion  coefficient  of  oxygen. 
The  initial  and  boundary  conditions  for  Eq.  (29)  are 

C?2 


/=o 


_ p^02 

“  S ,L 


(30) 


dC° 2 

epg(l  -  s)Dfe{{-f-  +  PgUgCf 


x=0 


— —  M°2 
4  F 


^y02  _  ^02 


L=  0 


(31) 


(32) 


2.2.5.  Heat  transfer 

Differentiating  Eqs.  (9)— (1 1),  rearranging  the  resultant  three 
equations  based  on  Eqs.  (1)  and  (2),  and  summating  them 
together,  the  final  energy  conservation  equation  becomes 

dT 

[sp\sc\  +  epg(l  -  s)cg  +  (1  -  s)pscs] 


dt 


+ 


9 

dx 


—(sk\  +  (1  -  s)kg  +  ks) 


dT 


dT 


dx 
h2o 


+\pmc\  +  PgWgCg]—  =  Hvapra  +  q 
where 

^vap  —  kg  —  h\ 


(33) 


(34) 


is  the  enthalpy  of  water  vaporization  and  generally  is  a  function 
of  temperature. 

The  external  volumetric  heat  or  sink  source  here  comes  from 
the  conduction  current  and  is  given  as 

2 


4  —  4\  +  4g  +  4s  —  ac 


(Jr 


(35) 


where  crc  is  the  electric  conductivity  of  solid  phase. 

Initial  condition  and  boundary  conditions  for  Eq.  (33)  are 


II 

O 

II 

(27) 

(sk  i  +  (1  - 

(28) 

ii 

II 

(36) 


dT 

dx 


=  lr  ~ 


T\x=qAS{ 


Pt 


x=0 


4  F 


0c 


(37) 


(38) 
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where  A  is  the  entropy  production  for  the  cathode  [13,14,21]. 
t]c  is  the  overpotential  in  the  cathode  catalyst  layer  and  related 
to  the  cathode  current  density  (zc)  by  [14] 

ic  =  £i((l  -  s)xg2 ) I *=0 exp( — ?2  fic )  (39) 

where  and  £2  are  two  parameters  adopted  from  [13,14],  x®2 
is  the  molar  fraction  of  oxygen. 

2.3.  Parameter  expressions  and  constitutive  relations 


D°2ff  =  £l5D°2 

g,eff  LJg 


(48) 


and  the  bulk  diffusivities  are  temperature  dependent  and  given 
by  [2] 


£>?2°  =  0.256  x  10"4 


T 


2.334 


X 


307.15 


(49) 


D?2  =  0.1775  x  10-4  x 


T 


1.823 


273.15 


(50) 


The  complexity  and  difficulty  of  the  numerical  analysis  for 
the  multiphase  transport  through  porous  GDL  arise  from  the 
fact  that  the  porous  GDL  properties  such  as  relative  perme¬ 
ability,  dynamic  viscosity  and  capillary  pressure  depend  upon 
liquid  water  saturation  and  temperature.  To  complete  the  needed 
numerical  study,  some  empirical  or  fitted  parameters  or  consti¬ 
tute  relations  are  employed  and  listed  here. 

Relative  permeabilities  of  liquid  phase  and  gas  phase  are 
related  to  the  liquid  water  saturation  by  [4], 

kr\  =  s3  (40) 

fcrg  =  (1  -  S)3  (41) 

The  dynamic  viscosity  of  liquid  phase  depends  on  temperature 
[14], 

IM  =  0.6612 (T  -  229)"1562  (42) 


A  fitted  relation  is  employed  to  reflect  the  dependence  of  the 
enthalpy  of  water  vaporization  //vap  (Jkg-1)  on  temperature 
based  on  the  data  provided  in  [22] . 

tfvap  =  3.0709  x  105(647.15  -  T)035549  (51) 

The  specific  heat  capacity  of  liquid  water  as  a  function  of  tem¬ 
perature  is  fitted  based  on  the  figure  provided  in  [23] 

d  =  -4.0699  X  10_8(T  —  273. 15)5  +  1.3113 

x  10“5(T  -  273. 15)4  -  1.6290  x  10_3(T  -  273.15)3 

+  1.0536  x  10_1(T  -  273. 15)2 

-3.2989(T- 273.15) +  4216.4  (52) 

The  specific  heat  capacities  for  oxygen,  nitrogen  and  water 
vapour  are  also  temperature  dependent  and  given  in  [23]  as 


The  capillary  pressure  expression  is  adopted  from  [10], 

pc  =  <rcos(0c)^  —  J  J(s)  (43) 

where  a  is  the  surface  tension,  0C  is  the  contact  or  wetting  angle, 
J(s)  is  the  Leverett  function  and  given  for  both  hydrophilic  and 
hydrophobic  porous  media  as  follows, 


— ^-(3.068  +  1.638  x  10“3T  -  0.512  x  10“6T2) 

M°  2 

(53) 

— *=-( 3.247  +  0.712  x  10_37’  -  0.041  x  10_6T2) 

Mn 2  v  3 

(54) 


1.417(1  -  s)  -  2.210(1  -  s)2  +  1.263(1  -  sf  if<9c  <  90° 
1.417s  -  2.210s2  +  1.263s3  if  <9C  >  90° 


(44) 


The  interfacial  mass  transfer  rate  of  water  between  the  gas  phase 
and  the  liquid  phase  is  given  by  [1,9,14,15] 

MH2° 

,H20  _  J  ^cond£(l  —  s')  — 


m  *  = 


Pl(YH20  _  H20 .  ifrH20  >  rH20 

vAg  Asat  )  11  Ag  —  Asat 


kvapsspipg(x^2°  -  x^2°) 


ifxj?2° 


<  V 


sat 

h2o 

sat 


(45) 


here  factors  (1  —  s)  and  s  are  introduced  to  account  for  the  varia¬ 
tion  of  the  interfacial  mass  transfer  with  liquid  water  saturation, 
kCond  and  kVap  are  the  water  vapour  condensation  rate  and  the 
liquid  water  evaporation  rate,  respectively,  x^ 2°  is  the  molar 
fraction  of  water  vapour,  x^2t°  is  the  saturation  molar  fraction 
of  water  vapour  and  related  to  the  temperature  by  [14], 

xh2o  (r  ^  atm)  =  io^28-59051_8-21ogio(r)+o.oo24804r-3i42.3i/r) 

(46) 


c"2°  =  — 77-7+3.634  +  1.195  X  10“3T  + 0.135  x  10“6T2) 

g  h2o  v  7 

(55) 

The  specific  heat  capacity  of  the  gas  phase  is  estimated  as 
cg  =  c°2C°2  +  Cg2°Cg2°  +  Cg2(l  -  C°2  -  Cg2°)  (56) 


The  effective  diffusion  coefficients  of  water  vapour  and  oxygen 
are  estimated  respectively  using  Bruggeman  relation, 


DU2°  =  £h5  DU2° 

g,eff  g 


(47) 


Similarly,  the  thermal  conductivity  of  the  gas  phase  is  estimated 
by 

r  —  *,02^02  ,  rH20^H20  ,  ilN2 /i  _  ^02  _  ^H20\ 

H  ~Kg^g  -h  Kg  cg  -h  Kg  G  cg  cg  1 


(57) 
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The  molar  fraction  of  each  species  is  related  to  its  mass  fraction 
by 


l  =  Cl 


;  Me 


Ml 


(58) 


E*g  = 1  <59> 

j 


where  j  —  O2,  H2O,  N2.  Therefore,  the  molar  mass  of  the  gas 
phase  can  be  calculated  by 

m*=(e§)  <6o» 

As  in  [14],  by  keeping  x^2/x^2  =21/79  at  the  cathode 
GDL/channel  interface,  the  mass  fractions  of  water  vapour  and 
oxygen  at  the  interface  can  be  determined  by  the  relative  humid- 


ity 

fx  =  RH  •  41° 

(61) 

1  -  rH2° 

xO,  _  1  Xg  ,L 

g’L  1+79/21 

(62) 

Therefore,  the  mass  fractions  of  water  vapour  and  oxygen,  and 
molar  mass  of  the  gas  phase  at  the  cathode  GDL/channel  inter¬ 
face  can  be  calculated  as 


rH  2° 
xg,L 


MH2° 


(63) 


(64) 


Mg,L  =  x°2t  M°2  +  xf2T°MU2°  +  (1  -  jeft  -  xf2T°)MN 2 


g  ,L- 


g  ,L  ,L 


(65) 


3.  Results  and  discussion 


The  parameter  values  used  in  the  base  case  are  listed 
in  Table  1.  In  order  to  investigate  the  effect  of  various 
parameters  on  the  transient  or  steady  behaviour  of  differ¬ 
ent  species  in  the  cathode  GDL,  four  parameters  (the  liq¬ 
uid  water  saturation  at  the  GDL/channel  interface  ( sl) ),  the 
percentage  of  liquid  water  to  all  the  water  at  the  CL/GDL 
interface  (/3),  the  current  density  (/c),  and  the  wetting  angle 
(0C))  are  chosen  to  be  varied.  The  time  needed  for  the  liquid 
water  saturation  to  reach  steady  state  is  defined  similar  as  in 
[24], 


/0L  s(t ,  x)  dx—  s(t+ At,  x)  dx 


•L 


^steady — max  sabs 


f0L  s(t ,  x)  dx 


<le— 5 


(66) 


where  At  is  the  time  step  used  in  the  numerical  simulation. 


Table  1 

Parameter  values  for  base  case 


e 

0.4 

P\  (kgm-3) 

983 

(kg  m-  1  s—  1 ) 

1.9  X  10“5  [14] 

K  (m2) 

7  x  10“13  [14] 

^cond  (s  1) 

100  [1] 

Lap  (N-1  m2  S-1) 

1/101325 

MH2°  (kg  mol-1) 

1.8  x  10-2 

M° 2  (kg  mol-1) 

3.2  x  10“2 

Mn 2  (kg  mol-1) 

2.8x  10“2 

R  (J  mol- 1  K-1) 

8.315 

so 

0 

SL 

0 

ic  (Am-2) 

10000 

F  (C  mol-1) 

96487 

a 

0.25  [14] 

pgL  (Nm-2) 

101325 

ch2° 

1 

0  (N  m- 1 ) 

6.25  x  10“2  [14] 

0C  (°) 

110  [14] 

RH  (%) 

95 

ki  (W  m- 1  K- 1 ) 

0.58  [25] 

k^1  (Wm-1  K"1) 

0.024  [25] 

(Wm-1  K-1) 

0.024  [25] 

^2°  (Wm-1  K-1) 

0.016  [25] 

ks  (Wm-1  K-1) 

1.7  [25] 

cs  (Jkg-'K-1) 

709  [25] 

ps  (kgm-3) 

2267  [25] 

Tl(  K) 

333 

AS*1  (J  mol- 1  K"1) 

-326.36  [14] 

crc  (Sm-1) 

1900  [14] 

?l  (Am"2) 

0.12  [14] 

?2  (V-1) 

30.6  [14] 

P 

0.5 

L  (m) 

3  x  i<r4 

3. 1 .  The  effect  of  sl 

The  parameter,  sl,  is  the  liquid  water  saturation  specified 
as  a  boundary  condition  at  the  GDL/channel  interface.  Fig.  2 
shows  the  change  of  the  distribution  of  the  liquid  water  satu¬ 
ration  in  the  GDL  with  time.  When  sl  =  0  (the  baseline  case), 
implying  no  liquid  water  exists  at  the  GDL/channel  interface, 
the  transient  behaviour  of  the  liquid  water  saturation  is  shown 
in  Fig.  2(a).  At  t-  0,  no  liquid  water  in  the  GDL,  that  is,  s{t , 
x)\t=o  =  0,  is  assumed.  Thereafter,  liquid  water  starts  to  accu¬ 
mulate  at  the  cathode  CL/GDL  interface  and  seeps  gradually 
from  the  CL/GDL  interface  (x  =  0)  to  the  GDL/channel  interface 
(x  =  L)  because  of  the  capillary  pressure.  It  takes  about  14.5  s  for 
liquid  water  to  reach  its  steady  state.  The  highest  water  satu¬ 
ration  is  about  0.077  at  the  cathode  CL/GDL  interface.  When 
*sx  =  0.05  and  5^=0  =  0,  liquid  water  comes  from  both  sides  of 
the  GDL  and  meets  at  a  point  close  to  the  GDL/channel  interface 
as  shown  in  Fig.  2(b).  Water  saturation  then  increases  across  the 
entire  GDL  thickness  and  reaches  its  steady  state  after  about 
1 1.59  s.  At  steady  state,  the  water  saturation  at  v  =  0  reaches  to 
0.081.  When  sl  is  specified  as  0.1  or  0.15  and  s\t=o  =  0,  similar 
to  the  case  of  sl  =  0.05,  liquid  water  seeps  first  from  both  sides 
of  the  GDL  and  then  meets  at  a  point  close  to  the  CL/GDL  inter¬ 
face  as  indicated  in  Fig.  2(c)  and  (d).  The  time  needed  for  the 
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Fig.  2.  The  effect  of  water  saturation  at  x  =  L  on  the  transient  behaviour  of  water  saturation  in  the  cathode  GDL.  (a)  sl  =  0,  (b)  sl  =  0.05,  (c)  sl  =  0. 1 ,  (d)  sl  -  0. 15. 


water  saturation  to  reach  steady  state  is  further  reduced  to  4.8 
and  2.1  s  and  the  water  saturation  at  x  =  0  is  0.109  and  0.154, 
respectively.  The  water  saturation  distribution  at  steady  state  is 
a  linear-like  function  through  the  thickness  of  the  GDL  and  a 
little  bit  higher  at  x  =  0  than  at  x-L  when  liquid  water  saturation 
is  higher  than  0.1,  which  is  shown  by  the  solid  line  in  Fig.  2(c) 
and  (d). 

By  comparing  Fig.  2(a)-(d),  it  can  be  seen  that  the  water 
saturation  at  the  cathode  CL/GDL  interface  increases  with  an 
increase  of  the  water  saturation  at  the  cathode  GDL/channel 
interface,  and  meanwhile,  the  time  needed  for  water  saturation 
to  reach  steady  state  decreases.  This  phenomenon  is  also  clearly 
reflected  in  Fig.  3. 

Fig.  4  gives  the  mass  transferrate  distribution  across  the  thick¬ 
ness  of  the  GDL  at  steady  state  for  the  four  cases,  sl  =  0,  0.05, 
0. 1  and  0.15.  mH2°  >  0  means  water  vapour  condensation  takes 
place  and  mH2°  <  0  implies  liquid  water  evaporation  occurs.  In 
the  baseline  case  (Fig.  4(a)),  water  condensation  takes  place  in 
the  most  portion  of  the  GDL  close  to  the  CL/GDL  interface 
side  due  to  high  water  vapour  concentration  in  this  region  (see 
Fig.  5(b)).  Water  evaporation  happens  only  in  a  small  portion  of 
the  GDL  close  to  the  GDL/channel  side  because  of  a  relatively 
low  water  vapour  concentration  within  this  side  (see  Fig.  5(b)). 
Fig.  4(b)-(d)  indicate  that  when  sl  value  goes  up,  the  mass  trans¬ 
fer  close  to  the  GDL/channel  interface  is  in  favour  of  evaporation 
and  the  rate  goes  down  accordingly. 

The  time  needed  for  water  vapour,  oxygen,  gas  phase  pres¬ 
sure,  and  temperature  to  reach  steady  state  are  usually  very  short 
(<0.2  s),  so  their  transient  behaviours  are  not  as  important  as  liq¬ 


uid  water  and  thus  omitted  here.  But  their  steady  state  behaviours 
are  presented  hereafter. 

Fig.  5  depicts  the  steady  state  distributions  of  gas  phase 
pressure,  water  vapour  concentration,  oxygen  concentration  and 
temperature  across  the  thickness  of  the  cathode  GDL  at  differ¬ 
ent  levels  of  liquid  water  saturation  at  x-L.  It  can  be  seen  that 
the  change  of  sl  value  has  a  slight  effect  on  the  gas  phase  pres¬ 
sure  drop  and  the  temperature  increase  within  the  GDL  from  the 
channel  side  to  the  CL  side,  as  shown  in  Fig.  5(a)  and  (d).  The 
gas  phase  pressure  drop  is  about  2.3  Pa  while  the  temperature 
increase  is  about  1.3  °C.  Fig.  5(b)  and  (c)  shows  that  an  increase 


s  (x=L) 


Fig.  3.  The  effect  of  water  saturation  at  x-L  on  the  time  needed  for  water 
saturation  to  reach  steady  state  and  on  the  water  saturation  at  x  =  0. 
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Fig.  5.  The  effect  of  water  saturation  at  x  =  L  on  the  steady  state  distribution  of  (a)  water  vapour  concentration,  (b)  oxygen  concentration,  (c)  gas  phase  pressure,  and 
(d)  temperature. 
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Fig.  6.  The  effect  of  parameter  ft  on  the  transient  behaviour  of  water  saturation  in  the  cathode  GDL.  (a)  ft  =  0,  (b)  f  =  0.25,  (c)  f  =  0.5,  (d)  f  =  1.0. 


of  sl  value  will  result  in  a  higher  water  vapour  concentration 
and  a  lower  oxygen  concentration  across  the  thickness  of  the 
cathode  GDL. 

3.2.  The  effect  of 

It  is  unclear  at  this  stage  that  how  much  liquid  water  and  how 
much  water  vapour  enters  the  GDL  from  the  cathode  CL  side. 
In  this  paper,  a  parameter  /3  is  introduced  to  account  for  the  liq¬ 
uid  water  percentage  to  all  of  the  water  at  the  cathode  CL/GDL 
interface.  Fig.  6  shows  the  effect  of  parameter  /3  on  the  transient 
behaviour  of  liquid  water  saturation  in  the  cathode  GDL.  If  all  of 
the  water  at  the  cathode  CL/GDL  interface  is  in  vapour  form,  that 
is,  /3  =  0,  the  water  vapour  concentration  is  greater  than  its  satu¬ 
ration  value  xfft°  and  as  a  result,  it  condenses  into  liquid  form. 
Therefore,  liquid  water  saturation  gradually  increases  from  the 
CL  side  to  the  channel  side  and  finally  reaches  its  steady  state 
where  liquid  water  pressure  and  gas  phase  pressure  are  balanced 
by  the  capillary  pressure  (pc  =  pg  —  pi).  The  whole  process  takes 
about  25.27  s  and  the  maximum  liquid  water  saturation  is  about 
0.058  at  the  cathode  CL/GDL  interface,  as  shown  in  Fig.  6(a). 
If  the  value  of  the  parameter  /3  is  respectively  increased  from 
0  to  0.25,  0.5  and  1.0,  meaning  that  more  liquid  water  enters 
the  GDL  from  the  CL  side,  liquid  water  saturations  at  steady 
state  are  also  increased,  as  shown  in  Fig.  6(b)-(d),  and  the  time 
needed  for  liquid  water  saturation  to  reach  steady  state  is  cor¬ 
respondingly  reduced  from  25.27  to  18.06,  14.56  and  10.76  s, 
respectively.  At  the  steady  state,  the  liquid  water  saturation  at  the 


CL/GDL  interface  associated  with  /3  =  0.25,  0.5  and  1.0  reaches 
0.069,  0.077  and  0.088,  respectively. 

The  variations  of  the  time  needed  for  liquid  water  saturation 
to  reach  steady  state  and  the  water  saturation  at  the  CL/GDL 
interface  with  the  parameter  (3  are  plotted  in  Fig.  7.  It  can  be 
seen  that  a  higher  /3  value  results  in  a  shorter  time  needed  for 
water  saturation  to  reach  steady  state  and  higher  liquid  water 
saturation  at  the  CL/GDL  interface. 

Fig.  8  gives  the  mass  transfer  rate  distributions  in  the  GDL  at 
steady  state  for  /3  =  0.0,  0.25,  0.5  and  1.0,  respectively.  Though 


Fig.  7.  The  effect  of  parameter  ft  on  the  time  needed  for  water  saturation  to 
reach  steady  state  and  on  the  water  saturation  at  x  =  0. 
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(a)  X  (m)  x  10" 


Fig.  8.  Mass  transfer  rate  at  steady  state,  (a)  f>  -  0,  (b)  f>  =  0.25,  (c)  ft  =  0.5,  (d)  f  =  1.0. 


the  interplay  among  phase  change,  pressure,  temperature  and 
boundary  conditions  is  very  complex,  some  quantitative  and 
qualitative  predictions  can  be  carried  out  using  the  set  of  models 
given  above.  Under  the  given  initial  and  boundary  conditions,  it 
can  be  seen  from  Fig.  8(a)-(c)  that  water  vapour  condensation 
usually  takes  place  near  the  CL  side  and  evaporation  occurs  near 
the  channel  side  because  of  the  water  vapour  concentration  is 
greater  than  its  saturation  value  near  the  CL  side  and  less  than 
its  saturation  value  near  the  channel  side.  The  dividing  point 
between  condensation  and  evaporation  moves  from  the  channel 
side  to  the  cathode  CL  side  when  (3  increases.  But  when  fi  is 
very  close  to  1,  say  0.95  </3<  1,  condensation  process  does  not 
happen  and  only  evaporation  takes  place,  as  shown  in  Fig.  8(d), 
because  too  much  liquid  water  exists  and  water  vapour  con¬ 
centration  through  the  whole  GDL  does  not  reach  its  saturation 
level. 

The  steady  state  distributions  for  gas  phase  pressure,  water 
vapour  concentration,  oxygen  concentration  and  temperature  in 
the  GDL  at  different  /3  values  are  given  in  Fig.  9.  Noting  that, 
when  /3  =  0  and  0.25,  the  gas  phase  pressure  at  the  CL/GDL 
interface  is  higher  than  the  pressure  at  the  GDL/channel  inter¬ 
face  due  to  the  higher  water  vapour  concentration  at  the  CL/GDL 
interface,  as  shown  in  Fig.  9(a).  The  gas  phase  pressure  begins 
to  drop  from  the  channel  side  to  the  CL  side  when  the  (3  value 
becomes  bigger  (/3  =  0.5  and  1  in  Fig.  9(a)).  The  water  vapour 
concentration  decreases  with  the  increase  of  fi  value,  as  shown 
in  Fig.  9(b),  which  leads  to  a  high  oxygen  concentration  distri¬ 
bution  (Fig.  9(c))  because  of  the  less  convection  effect  of  water 
vapour  on  the  oxygen  transport  in  the  GDL.  Furthermore,  it  can 


be  seen  from  Fig.  9(d)  that  a  bigger  /3  value  will  reduce  the 
temperature  difference  between  the  both  sides  of  the  GDL. 

3. 3.  The  effect  of  ic 

The  water  saturation  distributions  at  different  times  are  plot¬ 
ted  in  Fig.  10  for  four  current  densities,  ic  =  1000,  2000,  5000 
and  10,000  A  m~2.  A  higher  current  density  implies  more  water 
produced  through  the  electrochemical  reaction  in  the  cathode 
CL  and  more  oxygen  is  consumed.  Therefore,  more  liquid  water 
and  water  vapour  come  into  the  cathode  GDL  from  the  CL/GDL 
interface  that  results  in  a  high  liquid  water  saturation  distribution 
through  the  thickness  of  the  GDL,  as  shown  in  Fig.  10(a)-(d). 
The  maximum  liquid  water  saturation  for  each  current  density 
is  found  at  the  CL/GDL  interface  and  is  increased  from  0.039 
for  ic  =  1000  A  m-2  to  0.049,  0.064,  and  0.077,  respectively,  for 
ic  =  2000,  5000  and  10,000  Am-2.  The  time  needed  for  liquid 
water  saturation  to  reach  steady  state  is  down  from  76.37  s  for 
ic  =  1000  A  m~2  to  45.73,  23.42,  and  14.56  s  for  ic  =  2000,  5000 
and  10,000  Am-2,  respectively. 

The  effect  of  the  current  density  on  the  time  needed  for  water 
saturation  to  reach  steady  state  and  on  the  water  saturation  at 
the  CL/GDL  interface  are  shown  in  Fig.  1 1.  It  can  be  seen  that 
increasing  the  working  current  density  reduces  the  time  for  liq¬ 
uid  water  to  reach  steady  state  but  increases  the  liquid  water 
saturation  at  the  CL/GDL  interface. 

Fig.  12  gives  the  mass  transfer  rate  distributions  at  steady 
state  for  the  four  cases  corresponding  to  ic  =  1000,  2000,  5000 
and  10,000  A  m~2,  respectively.  It  can  be  seen  that  at  low  current 
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Fig.  9.  The  effect  of  parameter  f  on  the  steady  state  distribution  of  (a)  water  vapour,  (b)  oxygen,  (c)  gas  phase  pressure,  and  (d)  temperature. 


Fig.  10.  The  effect  of  the  current  density  on  the  transient  behaviour  of  water  saturation  in  the  cathode  GDL.  (a)  ic  =  1000Am  2,  (b)  ic  =  2000Am  2,  (c) 
ic  =5000  Am-2,  (d)  ic  =  10,000  Am-2. 
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Fig.  1 1 .  The  effect  of  the  current  density  on  the  time  needed  for  water  saturation 
to  reach  steady  state  and  on  the  water  saturation  at  x  =  0. 

density,  as  shown  in  Fig.  12(a)  and  (b),  evaporation  dominates 
the  phase  change  process  of  water  in  a  large  portion  of  the  GDL 
close  to  the  GDL/channel  interface  while  condensation  takes 
place  only  in  a  small  portion  of  the  GDL  close  to  the  CL/GDL 
interface.  The  reason  is  that  at  low  current  density  water  vapour 
enters  into  the  GDL  from  the  CL  side  and  is  oversaturated  in  a 
small  portion  of  the  GDL  close  to  the  CL/GDL  interface  side 
resulting  in  the  water  vapour  condensation  in  this  part.  On  the 
other  hand,  the  most  part  of  the  GDL  close  to  the  GDL/channel 
interface  side  is  unsaturated,  resulting  in  the  liquid  water  evap¬ 
oration.  With  the  increasing  of  the  current  density,  more  and 
more  liquid  water  enters  into  the  GDL  resulting  in  the  water 


vapour  in  the  most  part  of  the  GDL  from  the  left  side  to  the  right 
side  over  saturated,  and  gradually,  water  vapour  condensation 
replaces  evaporation  becoming  the  major  process  for  the  phase 
change  in  the  GDL,  as  shown  in  Fig.  10(c)  and  (d). 

The  distributions  for  gas  phase  pressure,  water  vapour  con¬ 
centration,  oxygen  concentration  and  temperature  through  the 
thickness  of  the  cathode  GDL  are  shown  in  Fig.  13.  It  can  be 
seen  that  increasing  the  current  density  will  increase  the  gas 
phase  pressure  drop  in  the  GDL,  and  lower  the  oxygen  concen¬ 
tration,  but  will  increase  the  water  vapour  concentration  and  the 
temperature  difference  between  the  two  sides  of  the  GDL. 

3.4.  The  effect  of  0 c 

Fig.  14  shows  the  dynamic  behaviours  of  liquid  water  sat¬ 
uration  in  the  GDL  versus  the  four  different  contact  angles, 
0C  =  0,  70,  110  and  150  °C.  The  first  two  cases  correspond  to 
hydrophilic  GDLs  and  the  last  two  represent  hydrophobic  GDLs. 
The  shape  of  the  water  saturation  distribution  across  the  thick¬ 
ness  of  the  GDL  seems  not  affected  too  much  by  the  contact 
angle  of  the  GDL  but  the  values  change  a  lot.  For  hydrophilic 
GDLs  (6>c  <  90  °C),  a  bigger  contact  angle  implies  that  the  GDL 
tends  to  hold  more  water  so  a  higher  steady  state  distribution 
of  liquid  water  saturation  is  expected,  Fig.  14(a)  and  (b);  for 
hydrophobic  GDLs  (0C  >  90  °C),  adverse  tendency  is  observed, 
i.e.,  increasing  the  contact  angle  will  reduce  the  ability  of  the 
GDL  to  hold  more  water  and  results  in  a  lower  steady  state  distri¬ 
bution  of  liquid  water  saturation,  Fig.  14(c)  and  (d).  In  addition, 
the  time  needed  for  water  saturation  to  reach  steady  state  will 


Fig.  12.  Mass  transfer  rate  at  steady  state  for  different  current  densities,  (a)  ic  =  1000  Am  2,  (b)  ic  =  2000  Am  2,  (c)  ic  =  5000  Am  2,  (d)  ic  =  10,000  Am  2. 
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Fig.  13.  The  effect  of  the  current  density  on  the  steady  state  distribution  of  (a)  water  vapour,  (b)  oxygen,  (c)  gas  phase  pressure,  and  (d)  temperature. 


Fig.  14.  The  effect  of  the  contact  angle  on  the  transient  behaviour  of  water  saturation  in  the  cathode  GDL.  (a)  6C  =  0  °C,  (b)  6C  =  70  °C,  (c)  6C  =  100  °C,  (d)  0C  =  150  °C. 
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Fig.  15.  The  effect  of  the  contact  angle  on  the  time  needed  for  water  saturation 
to  reach  steady  state  and  on  the  water  saturation  at  x  =  0. 

increase  when  the  contact  angle  increases  for  hydrophilic  GDLs 
while  the  time  will  decrease  when  the  contact  angle  increases 
for  hydrophobic  GDLs.  The  same  trend  is  also  observed  for  the 
liquid  water  saturation  at  the  cathode  CL/GDL  interface.  These 
phenomena  are  well  reflected  in  Fig.  15.  The  more  the  contact 
angle  closes  to  90  °C,  the  longer  the  time  needed  for  water  sat¬ 
uration  to  reach  steady  state  and  the  higher  the  water  saturation 
at  the  interface  of  the  cathode  CL/GDL  are.  The  steady  state 
mass  transfer  rate  and  distributions  of  water  vapour,  oxygen, 
gas  phase  pressure  and  temperature  are  almost  unchanged  for 
the  four  contact  angles  and  so  their  figures  are  not  presented 
here. 

4.  Conclusions 

The  transient  dynamics  of  water  transport  in  the  cathode 
gas  diffusion  layer  of  PEM  fuel  cells  is  studied  based  on  a 
one-dimensional,  two-phase,  nonisothermal  transient  model. 
According  to  the  numerical  simulations,  it  can  be  concluded 
as  follows: 

1 .  When  the  water  saturation  at  the  cathode  GDL/channel  inter¬ 
face  is  big  enough  (>0.1)  liquid  water  is  almost  uniformly 
distributed  across  the  thickness  of  the  cathode  GDL  at  steady 
state.  Increase  the  water  saturation  at  the  GDL/channel  inter¬ 
face  will  increase  the  entire  water  saturation  in  the  whole 
GDL  and  reduce  the  time  needed  for  liquid  water  saturation 
to  reach  steady  state.  In  order  to  obtain  a  high  oxygen  con¬ 
centration  at  the  CL/GDL  interface,  the  water  saturation  at 
the  GDL/channel  interface  should  be  kept  as  low  as  possible. 

2.  More  liquid  water  and  less  water  vapour  at  the  CL/GDL  inter¬ 
face  are  helpful  to  reduce  the  time  needed  for  liquid  water 


saturation  to  reach  steady  state  and  helpful  for  the  oxygen 
transport.  Water  vapour  condensation  takes  place  in  the  GDL 
near  the  CL  side  and  liquid  water  evaporation  occurs  near 
the  channel  side.  This  conclusion  depends  strongly  on  the 
boundary  conditions  at  the  GDL/channel  interface. 

3.  At  low  current  density  range,  evaporation  dominates  the 
water  phase  change;  but  at  high  current  density  range,  con¬ 
densation  is  the  major  process  for  water  phase  change  in  the 
cathode  GDL.  Increasing  the  current  density  will  reduce  the 
time  needed  for  water  saturation  to  reach  steady  state  and 
increase  the  water  saturation  at  the  CL/GDL  interface. 

4.  For  hydrophilic  GDLs  (0C  <  90  °C),  when  the  contact  angle 
increases  both  the  time  needed  for  liquid  water  saturation  to 
reach  steady  state  and  the  water  saturation  at  the  interface 
of  the  CL/GDL  increase  and  attain  their  maximum  point 
around  90  °C;  for  hydrophobic  GDLs  (0c>9O°C),  both  of 
them  decrease  when  the  contact  angle  increases. 
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